!// Dmitry.Konovalov@jcu.edu.au Dmitry.A.Konovalov@gmail.com, 2013

module func_mod; 
use memory_mod;
use grid_mod;
!use HYP_2F1_MODULE, only: HYP_2F1;
implicit none;
logical, save :: FUNC_TRACE = .FALSE.;

! FLF - func_log_fact
integer, save :: FLF_STEP = 100;
real(DP), pointer, save :: FLF_LOGS(:);


!????????
!type FUNC_T  !PT - pointer to type; or type-pointer
! real(DP), pointer :: x(:), y(:), dy(:), dy2(:);
! integer mini, maxi;
!end type
!
!type FUNC_VEC_T
!  real(DP), pointer :: x(:), y(:), drv(:), drv2(:); 
!end type
!type FUNC_VEC_PT  !PT - pointer to type; or type-pointer
!  type(FUNC_VEC_T), pointer :: fp; ! no better way 
!end type


contains 
subroutine test_func_mod(); 
  character(LEN=64) :: name = "test_func_mod";
  call dbg_in(name);
  call test_intrinsic ();
  call func_polynom_test();
  call test_func_fact();
  call test_cmplx_2f1();
  
  !ok = vec_dot_test_speed();
  call dbg_out(name); 
end subroutine


subroutine test_intrinsic(); 
  real(DP) :: maxerr;        
  character(LEN=64) :: name = "test_intrinsic";
  call dbg_in(name);
  
  call dbg_int("real(DP=", DP);
  call dbg_int("selected_real_kind(16)=", selected_real_kind(16)); 
  
  maxerr = EPS30;
  if (DP == 8) maxerr = EPS14;
  
  call dbg_assert_abs_diff_LE("EXP(0)", EXP(D0), D1, maxerr); 
  call dbg_assert_abs_diff_LE("EXP(1)", EXP(D1), EXP1, maxerr); 
  
  call dbg_assert_abs_diff_LE("LOG(EXP1)", LOG(EXP1), D1, maxerr); 
  
  call dbg_assert_abs_diff_LE("SIN(0)", SIN(D0), D0, maxerr); 
  !call dbg_assert_abs_diff_LE("calc_sin(PI)", calc_sin(PI), D0, EPS14); 
  call dbg_assert_abs_diff_LE("SIN(PI)", SIN(PI), D0, maxerr); 
  
  call dbg_assert_abs_diff_LE("COS(0)", COS(D0), D1, maxerr); 
  call dbg_assert_abs_diff_LE("COS(PI)", COS(PI), -D1, maxerr); 

  call dbg_out(name); 
end subroutine  

! done 18-Mar-2013
subroutine test_func_fact();         
  character(LEN=64) :: name = "test_func_fact";
  call dbg_in(name);
  call dbg_assert_abs_diff_LE("fact(0)", func_fact(0), D1, EPS13); 
  call dbg_assert_abs_diff_LE("log_fact(0)", func_log_fact(0), D0, EPS13); 
  call dbg_assert_abs_diff_LE("fact(1)", func_fact(1), D1, EPS13); 
  call dbg_assert_abs_diff_LE("fact(2)", func_fact(2), D2, EPS13); 
  call dbg_assert_abs_diff_LE("fact(3)", func_fact(3), D6, EPS13); 
  call dbg_assert_abs_diff_LE("fact(4)", func_fact(4), 24._DP, EPS13); 
  
  call dbg_assert_abs_diff_LE("fact(5)", func_fact(5), 120._DP, EPS13); 
  call dbg_assert_abs_diff_LE("log_fact(5)", func_log_fact(5), LOG(120._DP), EPS13);  

  call dbg_assert_abs_diff_LE("fact(6)", func_fact(6), 720._DP, EPS12); !!!NOTE ERR12 !!!!
  call dbg_assert_abs_diff_LE("log_fact(6)", func_log_fact(6), LOG(720._DP), EPS13);  

  call dbg_out(name); 
end subroutine  


!real(DP) function calc_sin(x) result(res);   
!  real(DP), intent(IN) :: x;   
!  real(16) dd;
!  dd = x;
!  res = SIN(dd);         
!!  complex(DP) c;
!!  c = CMPLX(D0, x);
!!  res = AIMAG(EXP(c));
!  return;
!end function


!! assert two functions are the same for the given X
subroutine func_assert_abs_diff_LE(info, x, f, f2, maxerr);  
  character(len=1024) :: mssg;
  character(LEN=*), intent(IN) :: info;
  real(DP), intent(IN) :: x(:), f(:), f2(:), maxerr;
  integer :: i, maxi, mini, min2, max2, maxf, minf;
  call vec_range(mini, maxi, x);
  call vec_range(minf, maxf, f);
  call vec_range(min2, max2, f2);  
  call dbg_assert_true("minf == mini", minf == mini);! the same start
  call dbg_assert_true("min2 == mini", min2 == mini);! the same start
  call dbg_assert_true("maxi <= maxf", maxi <= maxf);
  call dbg_assert_true("maxi <= max2", maxi <= max2);
  do i = mini, maxi
    write (mssg, *) "-->dbg_assert_abs_diff_LE_func(info=", info, ", i = ", i;
    call dbg_assert_abs_diff_LE(TRIM(mssg), f(i), f2(i), maxerr);
  end do
end subroutine

real(DP) function func_fact(n) result(res);   
  integer, intent(IN) :: n; 
  character(LEN=64) :: name = "func_fact";
  call dbg_in(name);          
  call dbg_assert_true("n >= 0", n >= 0);
  res = func_log_fact(n);
  res = EXP(res);
  return;
end function

real(DP) function func_log_fact(n) result(res);    
  integer, intent(IN) :: n; 
  character(LEN=64) :: name = "func_log_fact";         
  call dbg_in(name);
  call dbg_assert_true("n >= 0", n >= 0);
  if ((n+1) > SIZE(FLF_LOGS, 1))  call load_fact_logs(n + FLF_STEP);
  res = FLF_LOGS(n + 1);
  return;
end function

subroutine load_fact_logs(n);   
  integer, intent(IN) :: n; 
  integer i; 
  real(DP) di;    
  character(LEN=64) :: name = "load_fact_logs";         
  call dbg_in(name);
  call mem_del_vec(FLF_LOGS);
  FLF_LOGS => mem_vec(n+1);
  FLF_LOGS(1) = D0;
  do i = 1, n 
    di = real(i, DP);
    FLF_LOGS(i+1) = FLF_LOGS(i) + LOG(di);
  end do
end SUBROUTINE 

!-------------  func_polynom
  
subroutine func_polynom_test();  
  integer :: nv;
  real(DP) :: minv, maxv, res;
  real(DP), dimension(100) :: py;      
  real(DP), pointer :: px(:);
  real(DP), dimension(3) :: coeff;      
  character(LEN=64) :: name = "func_polynom_test";         
  call dbg_in(name);

  nv = 3;  minv = D01; maxv = 1.1D0;
  px => mem_vec(nv);  
  call grid_step(px, minv, maxv); 
  call vec_dbg("px=", px);
  call dbg_assert_abs_LE("3.px(1)=minv", px(1) - minv, EPS14);
  call dbg_assert_abs_LE("3.px(2)=minv", px(2) - (maxv+minv)/2., EPS14);
  call dbg_assert_abs_LE("3.px(nv)=maxv", px(nv) - maxv, EPS14);
  
  coeff = (/D1, D0, D0/);
  call fun_vec_polynom(py, px, coeff);  call vec_dbg("py=", py);
  res = D1;
  call dbg_assert_abs_diff_LE("1.py(1)", py(1), res, EPS14);
  call dbg_assert_abs_diff_LE("1.py(2)", py(2), res, EPS14);
  call dbg_assert_abs_diff_LE("1.py(3)", py(3), res, EPS14);

  coeff = (/D0, D05, D0/);
  call fun_vec_polynom(py, px, coeff); call vec_dbg("py=", py);
  res = D05;
  call dbg_assert_abs_diff_LE("05.py(1)", py(1), res * px(1), EPS14);
  call dbg_assert_abs_diff_LE("05.py(2)", py(2), res * px(2), EPS14);
  call dbg_assert_abs_diff_LE("05.py(3)", py(3), res * px(3), EPS14);

  coeff = (/-D05, -D1, D3/);
  call fun_vec_polynom(py, px, coeff); call vec_dbg("py=", py);
  res = -D05;
  call dbg_assert_abs_diff_LE("513(1)", py(1),res-px(1)+3*px(1)*px(1),EPS14);
  call dbg_assert_abs_diff_LE("513(2)", py(2),res-px(2)+3*px(2)*px(2),EPS14);
  call dbg_assert_abs_diff_LE("513(3)", py(3),res-px(3)+3*px(3)*px(3),EPS14);
  
  call mem_del_vec(px);
end subroutine


!
SUBROUTINE fun_vec_polynom(res, x, coeff);     
  real(DP), intent(IN)  :: coeff(:), x(:);      
  real(DP), intent(OUT) :: res(:);      ! vector-x, vector-y
  integer :: i, mini, maxi;
!  character(LEN=64) :: name = "fun_vec_polynom";         
!  call dbg_in(name);
  call vec_init_res(res, x);
  call vec_range(mini, maxi, x); 
  do i = mini, maxi 
    res(i) = func_polynom(x(i), coeff);
  end do
end SUBROUTINE

SUBROUTINE vec_sin(res, x);     
  real(DP), intent(IN)  :: x(:);      
  real(DP), intent(OUT) :: res(:);      ! vector-x, vector-y
  integer :: i, mini, maxi;
  call dbg_str("-->vec_sin()", FUNC_TRACE);
  call vec_init_res(res, x);
  call vec_range(mini, maxi, x); 
  do i = mini, maxi 
    res(i) = SIN(x(i));
  end do
end SUBROUTINE

SUBROUTINE vec_cos(res, x);   
  real(DP), intent(IN)  :: x(:);      
  real(DP), intent(OUT) :: res(:);      ! vector-x, vector-y
  integer :: i, mini, maxi;
  call vec_init_res(res, x);
  call vec_range(mini, maxi, x); 
  do i = mini, maxi 
    res(i) = COS(x(i));
  end do
end SUBROUTINE

SUBROUTINE vec_exp(res, x);    
  real(DP), intent(IN)  :: x(:);      
  real(DP), intent(OUT) :: res(:);      ! vector-x, vector-y
  integer :: i, mini, maxi;
  call vec_init_res(res, x);
  call vec_range(mini, maxi, x); 
  do i = mini, maxi 
    res(i) = EXP(x(i));
  end do
end SUBROUTINE

real(DP) function func_polynom(x, coeff) result(y);   
  real(DP), dimension(:), intent(IN) :: coeff;      
  real(DP), intent(IN) :: x;      
  real(DP) :: powx;      
  integer :: i, mini, maxi;
!  character(LEN=64) :: name = "func_polynom";         
!  call dbg_in(name);
  call dbg_assert_true("SIZE(coeff, 1) > 0", SIZE(coeff, 1) > 0);
  call vec_range(mini, maxi, coeff);
  powx = D1;
  y = coeff(mini);
  do i = mini + 1, maxi
    powx = powx * x;
    y = y + powx * coeff(i);
  end do
  return;
end function


subroutine test_cmplx_2f1();         
  complex(DP) a, b, c, z, res, corr;
  complex(DP) :: HYP_2F1, cdp;
  character(LEN=64) :: name = "test_cmplx_2f1";
  call dbg_in(name);   
  
  a = ZEROc;  b = ZEROc;  c = ONEc;  z = ZEROc;
  res = func_cmplx_2f1_SLOW(a, b, c, z);  
  call dbg_assert_cabs_diff_LE("2F1(0,0,0,0)", res, ONEc, EPS30); 
  
!    // F_21(a, b; c; z) = F_21(1, 1; 2; z) = -z^(-1) ln(1-z)   // p.370 Abr&Steg
!    Cmplx a = new Cmplx(1, 0);
!    Cmplx b = new Cmplx(1, 0);
!    Cmplx c = new Cmplx(2, 0);
!    double zRe = 0.5;
!    Cmplx z = new Cmplx(zRe, 0);
!    res = Cmplx2F1.calc(a, b, c, z, Calc.EPS_10);
!    double corr = -Math.log(1. - zRe) / zRe;
!    log.assertZero("F_21(1, 1; 2; z=0.5)=", corr - res.getRe(), 2e-10);
  a = ONEc; b = ONEc; c = TWOc; z = CMPLX(D05, D0);  
  res = func_cmplx_2f1_SLOW(a, b, c, z);  
  corr = -LOG(ONEc - z) / z;  
  call dbg_assert_cabs_diff_LE("2F1(1,1,2,0.5)", res, corr, EPS15); 

  res = HYP_2F1(a, b, c, z);  
  call dbg_assert_cabs_diff_LE("HYP_2F1(1,1,2,0.5)", res, corr, EPS15); 
  
  !//http://mathworld.wolfram.com/HypergeometricFunction.html
  !// 2F1(1,1;1;z) = 1/(1-z)
  a = ONEc; b = ONEc; c = ONEc; z = CMPLX(D05, D0);    
  res = func_cmplx_2f1_SLOW(a, b, c, z);  
  corr = ONEc / (ONEc - z);  
  call dbg_assert_cabs_diff_LE("2F1(1,1,1,0.5)", res, corr, EPS30); 
  res = HYP_2F1(a, b, c, z);  
  call dbg_assert_cabs_diff_LE("HYP_2F1(1,1,1,0.5)", res, corr, EPS15); 
  
  a = ONEc; b = ONEc; c = ONEc; z = CMPLX(D0, D05);    
  corr = ONEc / (ONEc - z);   call dbg_cmplx("corr=", corr);
  res = func_cmplx_2f1_SLOW(a, b, c, z);  call dbg_cmplx("res =", res);
  call dbg_assert_cabs_diff_LE("2F1(1,1,1,0.5i)", res, corr, EPS30); 
  res = HYP_2F1(a, b, c, z);   call dbg_cmplx("res =", res);
  call dbg_assert_cabs_diff_LE("HYP_2F1(1,1,1,0.5i)", res, corr, EPS15); 
  
end subroutine  


complex(DP) function func_cmplx_2f1_SLOW(a, b, c, z) result(res);   
  complex(DP), intent(IN) :: a, b, c, z;      
  complex(DP) :: a2, b2, c2, n2, newRes, t;
  character(LEN=64) :: name = "func_cmplx_2f1_SLOW";
  call dbg_in(name);   
  call dbg_assert_true("ABS(z)<1", ABS(z) < D1);
  call dbg_assert_true("c.NOT.neg or zero int", .NOT. dbg_neg_int_or_zero(c));
  n2 = ONEc;
  res = ZEROc;  
  newRes = ONEc;  
  a2 = a;
  b2 = b;
  c2 = c;
  t = ONEc;
  do while (newRes .NE. res)  !// STEP-1: newRes = 1; res = 0; n = 1
    t = t * a2 * b2 * z / (c2 * n2) !// 1:t=a*b/c *z/1; 2:t=a*b/c *z/1  * (a+1)*(b+1)/(c+1) *z/2
    res = newRes; !// STEP-1: res = 1;
    a2 = a2 + ONEc; !// STEP-1: a2=a+1;
    b2 = b2 + ONEc; !// STEP-1: b2=b+1;
    c2 = c2 + ONEc; !// STEP-1: c2=c+1;   
    n2 = n2 + ONEc; !// STEP-1: n2=n+1;
    newRes = res + t;
    !call dbg_cmplx("res=", res);
    !call dbg_cmplx("new=", newRes);
  enddo
  call dbg_out(name);
end function


end module 